Basin record of a Miocene lithosphere drip beneath the Colorado Plateau

The sinking of gravitationally unstable lithosphere beneath high-elevation plateaus is proposed to be a key driver of their uplift. Numerical geodynamic models predict that lithosphere removal can lead to transient, dynamic topographic changes that could be preserved in the surface record, particularly in sedimentary deposits of lakes or playas that are subsequently inverted. However, few such examples have been documented. Here we show that the Miocene Bidahochi Basin, which was partially and intermittently filled by the Hopi Paleolake, preserves a record of the quasi-elliptical surface response to a viscous drip of lithosphere >100 km beneath the Colorado Plateau. New detrital zircon U-Pb, Lu-Hf, and trace-element data reveal systematic isotopic, geochemical, temperature and fO2 transitions in magmatism proximal to the basin. Integration of geophysical, geochemical, and geological evidence supports a spatially and temporally varying record of subsidence and uplift that is consistent with models of progressive dripping beneath plateaus with thick lithosphere. We demonstrate that dynamic topography at the scale of individual lithosphere drips can be recognized on the Colorado Plateau, despite the strength of its lithosphere.

asymmetry. At the same time, lithospheric drip(s) (with negligible horizontal component of advection) begin to form further from the edge of the plateau. Our model for the Miocene lithosphere drip beneath the Bidahochi Basin focuses at the scale of these individual drip and therefore does not account for the lithospheric step at the larger, regional scale. The Bidahochi Basin was likely further inboard from the Miocene edge of the plateau, not directly above a lithospheric step. The surface evidence in the Bidahochi Basin also does not point to an asymmetrical, laterally (towards the center of the plateau) migrating basin. For these reasons we do not further consider this model here.

fH2O
The design of the numerical model used in this paper also incorporates wet flow laws to take into account the likelihood that the upper mantle is partly hydrated, based on evidence from fourier transform infrared spectroscopy of xenoliths from the base of the Colorado Plateau mantle lithosphere as well as evidence that the upper mantle likely has up to 200 ppm of water (~45 ppm in olivine, Li et al., 2008;Hirschmann, 2006). 13,14 Uncertainties in the interpretation of geological, geophysical, and geochronologic data and geodynamic models

Geochronology
There are several sources of uncertainty in the use of isotope geochemistry and trace and rare earth elements as proxies for magmatic evolution. An inherent problem for detrital zircon geochronology is that the precise magmatic composition from which the zircons crystallized is unknown, and the partitioning behavior of various elements remains the subject of ongoing research, due to complexities in crystal interface/melt kinetics, the effect of fO2, and possible disequilibrium behavior (Claiborne, 2006) 15 . Nevertheless, the use of isotopic tracers (Hf) and other proxies like Ti concentration and U/Yb for fractionation or as indicators of melt source is well established. Use of Ti-in-zircon thermometer in detrital zircons is dependent on unknown aSiO2 and aTiO2 (Ferry and Watson, 2007) 16 . For this reason, we show an extended error bar that encompasses a large range of aTiO2 (Watson and Harrison, 2005) 17 , and further note that the effect of subunity aSiO2 and aTiO2 compete with each other. For example, in the case of aSiO2 = 0.5 and aTiO2 = 0.5, there would be no difference in the estimated crystallization temperature. Because the melts from which the zircons originated were likely silica undersaturated, aSiO2 is most likely <1. These uncertainties preclude the determination of an absolute crystallization temperature for any given grain; nevertheless, the trend of Ti concentration data reflects a systematic evolution of the melt source.
The calibrated oxybarometer of Loucks et al. 2020 18 uses ratios of U, Ce, and Ti without the independent determination of their ionic charge, the crystallization temperature, and melt composition. The oxybarometer is derived from thermodynamic relations and empirically calibrated with known zirconwhole rock pairs (Loucks et al., 2020) 18 . An associated uncertainty of 0.6 log unit fO2 is shown in Fig. 3. Additional uncertainty arises from the large analytical uncertainties associated with measurements of LREE (particularly La and Pr, which are <1ppm in many of the zircon analyses, as well as Nd and Sm), so we have avoided reliance on those LREE with the largest standard deviation (based on analyses of external standards with known REE concentrations).
We consider the trends shown in Fig. 3 to be robust despite the uncertainties discussed above, as the trends are supported by the correlation of each of these datasets: U-Pb date, fO2 proxy, inferred crystallization temperature, Hf isotopic signature, Hf/Lu, and U/Yb. The use of Hf/Lu ratio as a fractionation proxy, however, is less well established, because zircon Hf/Lu does not necessarily reflect the Hf/Lu ratio of the melt. Progressive zircon crystallization depletes Hf with respect to Lu in the melt, because Hf is preferentially incorporated in zircon (on the order of 10,000 ppm). However, the Hf concentration in zircon can be controlled by Zr/Hf in the melt (Grimes et al., 2015) 19 . Zr/Hf is relatively constant in most melts (~35-45, close to the chondritic ratio) because of their nearly identical behavior due to their 4+ charge and close ionic radii (Claiborne et al., 2006) 15 , but studies have demonstrated that in felsic (and particularly leucogranitic) melts, zircon Hf composition increases with greater magmatic fractionation, despite the decrease in Hf in the melt, because Zr/Hf can decrease to as low as 15.
Additional uncertainty (beyond analytical uncertainty) is associated with the conversion of trace element concentration in zircon to estimated whole rock values using partition coefficients. For this reason, none of the datasets shown involve calculations using partition coefficients, except  Fig. S6 shows the uncertainty for each analysis, propagated from both alpha and beta values for each partition coefficient as well as the greatest standard deviations of the external standards analyzed with the unknowns. The inset in Fig. S6 shows two examples of known whole-rock-zircon pairings from Chapman et al. (2016) as a reference for the expected scatter in zircon REE analyses using the same partition coefficients and analytical methods. The large uncertainties do not permit unambiguous identification of a source, but suggest that the REE geochemistry of the zircons are consistent with derivation from either source.

Geological
Sourcing of the youngest Bidahochi detrital zircon population A source of uncertainty in the attribution of young (~9-6 Ma) detrital zircons to proximal magmatic centers is that detrital zircons do not provide definitive information on their source. However, there are limited late Miocene magmatic centers that can plausibly have been the source of these zircons. Besides the HBVF, these include the olivine phyric hawaiites, trachyandesites, and trachytes of the Mt. Baldy volcanics (Nealey, 1988) 21 and the ~7-4 Ma isotopically evolved minette dikes (Placerville and Lizardhead dike swarms)(Lake and Farmer, 2015) and 10-9 Ma felsic porphyry of the western San Juans (WSJ, Gonzales, 2015) 22 . Given the scarcity of 9-6 Ma zircons on the Colorado Plateau and Colorado River detrital record (<0.1%) in general, it would be extremely unlikely for zircons of that age to be a significant population (~2% in the Bidahochi Fm) unless the source was proximal 23 : Precisely because zircons are modally minor in silica-understurated and alkaline volcanic rocks, the probability of finding them far from their source becomes vanishingly small outside an endorheic basin, because any mixing with other felsic, volumetrically-dominant sources on the Colorado Plateau with orders-of-magnitude greater concentration of zircon of other ages would drastically decrease the probability of encountering them in strata outside an internally-drained basin.
Additionally, though zircon crystallization is uncommon in silica-undersaturated volcanic rocks, there is ample evidence of zircon occurrence in basalts, including silica-deficient ocean island basalts, which may co-crystallize coeval badellyite and zircon (Grimes et 27 . These Zr concentrations correspond to around 0.1 wt% normative zircon. Though zircon grains may be modally minor in such rocks, widespread intrabasinal dispersion and redistribution of volcaniclastic detritus from pyrcolastic fallout, vessicular tephra, bombs, and volcanic flows would easily concentrate them in an endorheic basin after the weathering of large volumes of volcanic material over geologic time.
Third, the population of 9-6 Ma zircons from the Bidahochi Formation in this study are more consistent with the K/Ar geochronology and other data from the HBVF and Mt. Baldy volcanic rocks than dates from the more distal sources such as the San Juans (Gonzales, 2015) 22 (Supplementary Fig. 1). In particular, sourcing from the 10-9 Ma felsic or isotopically evolved rocks of the San Juans is inconsistent with the isotopic composition of the ~9-8 Ma Ma zircons. On the other hand, the Bidahochi detrital age spectrum mirrors the K/Ar dates reported in Dallegge (1999) 28 , including a subtle (at the Ma-resolution) bimodal distribution of dates around ~8 Ma and 6.5-7 Ma (Supplementary Fig. 1). Notably, nearby deposits on the Black Mesa that are potentially correlative to the upper Bidahochi Fm, and that are almost certainly derived from the San Juan Volcanic Field volcaniclastic apron (based on high resolution K/Caage correlations of discrete 29-26 Ma ignimibrite events and detrital sanidine populations) ( 20 . Given the analytical uncertainties and that of the partition coefficients discussed above, the composition of the young zircons are consistent with derivation from either HBVF or Mt. Baldy, but cannot distinguish between these two sources.

Geometry of the Bidahochi Basin
The uncertainty regarding whether the zircon in the Bidahochi Basin were sourced from the Mt. Baldy Volcanics or the HBVF is relevant for our estimates of the size of the basin, and the wavelength of the sublithospheric instability that would have been responsible for the subsidence. If the HBVF and Mt. Baldy volcanics, both of which are alkalic and similarly high in Zn/Fe, were melts generated from flanking upwelling loci on either side of the basin, it would imply a larger basin (c. 250-300 km) than suggested by the extant outcrops of the lower Bidahochi Fm (c. 150 km). Such a larger basin is supported by the possible correlation of the Fence Lake Formation in New Mexico to the upper Bidahochi Fm. This range (from 250-300 km) corresponds to the upper end of the range of wavelength considered in our analysis in Fig. 5b.
The shape of the basin, as defined by the unconformable contact between the lower member and underlying Mesozoic rocks, has been considered by several authors (Dickinson et al., 2013;Karlstrom et al, 2017). 29,30 Various interpolation techniques and datasets have led to slightly different basin shapes (Fig. S7), but a commonality in all models is that the unconformity beneath the thickest lacustrine deposits (further to the northwest) sits at a higher elevation than the same unconformity to the southeast. The cross section in Fig. 4c is constrained by the control points listed in Supplementary Table 1. Fig. S8 shows an annotated version of the cross section, with superimposed stratigraphic constraints, and alternative interpretations of the onlap relationship between the Bidahochi strata and the basin margin that are permitted by these constraints. At the scale relevant to this study, there is no discernable difference between these interpretations.

Free air anomaly
The use of free-air anomaly to discern dynamic topography (i.e. topography not in isostatic equilibrium) has several caveats (Molnar et al., 2015) 31 . The first is that even isostatically compensated topography has free-air anomaly. Isostasy assumes the balancing of excess or deficient mass at the surface by mass at depth r, but the contribution to gravity is proportional to 1/r 2 . Thus, the possibility that some density anomaly (e.g. the slightly less dense Mesozoic rocks in the Kaiporowitz/Henry Mountain basins) in the lithosphere contributes to the gravity anomaly cannot be excluded. However, the free-air anomaly due to any isostatically compensated density anomaly of a geologically plausible magnitude would be near-zero (Chase et al., 2002) 32 . The magnitude of free-air anomalies further decreases as compensation depth decreases (relative to the relevant wavelength of topography/density anomaly) (Molnar et al., 2015) 31 . There is no indication of any correlation between the thickness of extant Mesozoic and younger strata and free-air anomaly, suggesting that the mass deficit from slightly-less-dense Mesozoic strata is largely compensated (Supplementary Fig. 4). Even if the density anomalies are completely uncompensated, they would contribute on the order of only ~4-7 mgal/km of Mesozoic strata, assuming a near-surface density deficit of ~80-170 kgm -3 (Ander 1981, Langenheim et al., 2000 33,34 . The second caveat is that the 138 mGal/km conversion assumes, conservatively, that the contribution to gravity of the density variation between the downwelling lithosphere and asthenosphere is negligible (because of the 1/r 2 relationship). To the extent that this assumption is not true, the potential dynamic subsidence implied would be even greater. Finally, free-air anomaly can result from flexural support of short-wavelength topography. The free-air gravity anomaly that is attributable to flexurally supported short-wavelength topography differs depending on the flexural-isostatic model used. In this case, we directly used the GRACE satellite missions data available in GeoMapApp, instead of applying a low-pass filter (Supplementary Fig. 3). This can also be compared to other studies that have explicitly considered the flexural-isostatic support of topography and/or the flexural-isostatic response to incision on a regional scale using a variable elastic thickness model (cf. Lazear

Tomography
The tomographic anomaly beneath the Escalante region and others elsewhere in the Colorado could be or have been variously interpreted as "lithospheric dripping," "downwelling," "foundering," "convective removal," "delamination," or "delamination-style convective lithospheric downwelling." We use the term lithospheric drip as a descriptive, non-genetic term to describe the apparent shape of the velocity anomaly, though the resolution of tomography may obscure the actual size and shape of the feature. For this reason we draw a series of contours to portray the seismically fast feature (Fig. S10). The term "lithospheric drip" emphasizes the viscous rheology that is dominant at that temperature and pressure range, and the small-scale nature of the anomaly as opposed to the wholesale removal of lithosphere, but does not necessarily imply the dominance of any particular process (e.g. convective removal, negative buoyancy, melt invasion).

Viscosity and timescale
The sequence of the initial increase in shear stress, prolonged maintenance of elevated shear stress, and final acceleration in advection velocity is a commonality in all our drip models, though the timescale of the process differs depending on a number of under-constrained parameters of the lithosphere and asthenosphere, including the mantle potential temperature, water fugacity, and density contrast. For example, a decrease in just 20 kgm -3 in density contrast, holding all other parameters constant, more than doubles the timescale the instability formation (Fig. S9). Note that the duration from initial onset of subsidence to the final sinking of the drip shown in Fig. S9 can also be in part attributed to the progression of adjacent drips, which also contribute to the stress distribution in the lithosphere. In a series of experiments, we also varied the water fugacity of the mantle lithosphere to control its viscosity, keeping the fH2O of the mantle asthenosphere constant (Fig. S5). Decreased viscosity or increased density contrast shortens the timescale of instability formation. Numerous other papers on numerical modelling have extensively tested and discussed at length the effect of these varying parameters (Gogus and Psyklywec, 2009; Gogus et al., 2022; Gogus et al., 2017), and we refer interested readers to those contributions. Because it is easy to significantly change the timescale of the models by adjusting any of these parameters (Fig. S9), we do not consider the absolute timescale in our interpretation of these models and instead focus on relative changes that are common throughout all models.